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ABSTRACT 



We investigate the formation of spatial structure in dense, self-gravitating particle systems 
such as Saturn's B-ring through local A'"-body simulations to clarify the intrinsic physics based 
on individual particle motion. In such a system, Salo (1995) showed that the formation of 
spatial structure such as wake-like structure and particle grouping (clump) arises spontaneously 
due to gravitational instability and the radial velocity dispersion increases as the formation 
of the wake structure. However, intrinsic physics of the phenomena has not been clarified. 
We performed local iV-body simulations including mutual gravitational forces between ring 
particles as well as direct (inelastic) coUisions with identical (up to A?" ~ 40000) particles. In 
the wake structure particles no longer move randomly but coherently. We found that particle 
motion was similar to Keplerian motion even in the wake structure and that the coherent 
motion was produced since the particles in a clump had similar eccentricity and longitude 
of perihelion. This coherent motion causes the increase and oscillation in the radial velocity 
dispersion. The mean velocity dispersion is rather larger in a more dissipative case with a 
smaller restitution coefficient and/or a larger surface density since the coherence is stronger in 
the more dissipative case. Our simulations showed that the wavelength of the wake structure 
was approximately given by the longest wavelength = 47r^G'S/K^ in the linear theory of 
axisymmetric gravitational instability in a thin disk, where G, E, and k are the gravitational 
constant, surface density, and a epicyclic frequency. 
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1 Introduction 



The observation by the Voyager revealed that Saturn's ring is not homogeneous but it has 
complex structures. The ring has fine axisymmetric subrings with the width ~ 10 km (e.g., 
Smith et al. 1982). Furthermore, the observation of Photopolarimeter (PPS) stellar occultation 
showed the finer density structure (microstructure) with scale down to 100m in B-ring (Espos- 
ito et al. 1983). Several ideas have been proposed to explain why such structure is formed 
and maintained against viscous diffusion due to collisions and gravitational scatterings of ring 
particles, but there is no satisfactory explanation. 

To account for the axisymmetric subring structure, viscous instability was proposed by Lin 
and Bordenheimer (1981), Ward (1981), and Lukkari (1981). If viscous instability arises, dif- 
fusion from a less dense region is larger than that from a denser region so that the density 
contrast is strengthen until non-linear effect becomes important. The condition of this instabil- 
ity is 9 (i/E)/9E < 0, where v and E are kinematic viscosity and surface density of ring. It has 
been investigated whether a ring particle system satisfies the condition or not. However, neg- 
ative results for viscous instability were reported by both theoretical studies using Boltzmann 
equation (Araki and Tremaine 1986, Araki 1988, 1991) and local A^-body simulations (Wisdom 
and Tremaine 1988, Salo 1991, 1992a, Richardson 1994). 

On the other hand, Salo (1992b), Richardson (1994), and Salo (1995) gave important sug- 
gestions on the microstructure observed in B-ring. They extended the local A-body simulations 
developed by Wisdom and Tremaine (1988) to include self-gravitational forces between parti- 
cles as well as inelastic coUisions and used larger number of particles (A^ > 1000) than that 
in Wisdom and Tremaine (1988) (A" = 40). They found the wake-like structure which may 
correspond to the substructure. Salo (1995) showed that the wakes arc created if both self- 
gravity and inelastic collisions are included. Only self-gravity or only inelastic collision does 
not create such structure. The wake structure is time dependent and transient, being created 
and destroyed on the time scale of an order of Keplerian period. The wake structure looks like 
fluid turbulence. 
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Salo (1995) also found that when the wake-hke structure is created, the equihbrium veloc- 
ity dispersion of ring particles always satisfies Q ~ 2, where Q is Toomre's non-dimensional 
parameter (Toomre 1964) defined by 

s.secE' 

where G is the gravitational constant, k is a epicyclic frequency, and is radial velocity 
dispersion. 

Since the wake structure does not develop when Q > 2, self-gravitational instability would 
be responsible for the wake structure. The linear calculation shows that a self-gravitating, 
differentially rotating disk of coUisionless particles becomes unstable against axisymmetric per- 
turbation for Q < 1 (e.g., Toomre 1964, Julian and Toomre 1966). For non- axisymmetric 
perturbation, Griv (1998) derived a similar criterion, Q ~ 2^1q/k, where is the angular 
velocity. Griv (1998) also confirmed the above criterion by local A^-body simulations. 

Since the relation of Q ~ 2 is maintained once the wake-like structure is formed (Salo 1995), 
the radial velocity dispersion increases with E. This increase in the radial velocity dispersion 
must be closely related with the formation of the wake structure. Salo (1995) suggested that 
such an increase of the velocity dispersion would come from scatterings by the collective wakes. 
He also suggested that the systematic motion of particles in the wakes may be responsible for 
the increase in the radial velocity dispersion. However, it has not been clarified how the radial 
velocity increase is related with the formation of the wake structure. 

The study on the wake structure is important to clarify the microstructure in Saturn's B- 
ring. Furthermore, it is also important in the issue of the stability of Saturn's B-ring, since 
the transient wake structure would induce relatively large angular momentum flux to quickly 
diffuse out the ring. (The detailed analysis of the angular momentum transfer will be done in 
the next paper). 

We perform local iV-body simulations including both mutual gravitational forces between 
particles and inelastic coUisions as Salo (1995) did. We focus on the problem of the relation 
between the wake formation and the velocity increase through detailed analysis about how 
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self-gravitating particles behave in the wake-like structure. 

Simulation method is described in section 2. In section 3, we first compare our results with 
previous works (especially with Salo 1995), then we will analyze the coherent particle motion in 
more detail. We find regular oscillation of the velocity dispersion associated with the increase 
in the time averaged velocity dispersion and clarify the relationship between the wake structure 
and the velocity change. We also analyze the scale of the structure by performing simulations 
with up to ~ 40000 particles. In section 4 we summarize the results and give discussion. 
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2 Numerical methods 
2.1 Model description 

We adopt local iV-body method which was first applied for the study of a dense ring system by 
Wisdom and Tremaine (1988) and followed by Salo (1991, 1992a, 1992b, 1995) and Richardson 
(1994). The "local" means that we consider a box with width and height Ly at a semimajor 
axis ao in the ring, which revolves in a circular orbit with the Kepler angular velocity VLq at the 
reference point oq and is small compared to the width of the whole ring (see Fig. 1). Motion 
of particles is pursued only in this box with the periodic boundary conditions. This method 
would be valid because we are considering the structure with much smaller scale than the width 
of the B-ring and the orbits of ring particles are nearly circular in the coplanar plane. 

We use the Cartesian coordinates at ao which rotate with the angular velocity VLq, in which 
X-axis points radially outward, y-axis points to the direction of orbital motion of the coordinate 
origin, and 2;-axis points to the direction normal to the orbital plane. Motion of a ring particle 
i is described by Hill's equation (e.g.. Hill 1878, Nakazawa and Ida 1988): 

(2) 

Zi), 

where irij and rjj are mass of particle j and the relative distance between particle i and j, 
and VLq = ^GMjal, Ms is mass of a central planet. The first and second terms on the 
right hand sides of Eqs. (0) denote the Coriolis and the tidal force, and the last terms are the 
mutual gravitational force between ring particles. When the mutual gravitational force can be 
neglected, Eqs. (0) are solved analytically and the solution represents an orbit of the Keplerian 
motion (see Nakazawa and Ida 1988). 

Scaling and parameters which characterize a ring system are as follow. Throughout the 
present paper, we assume that all particles have the same radius rp and mass m. The mutual 
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Hill's (tidal) radius is defined by rh = ha^, where h is the reduced Hill's radius defined as 



(3) 



rh represents the radius of potential wall between two particles at which central force (tidal 
force) and mutual gravitational force between the particles are balanced. Within the radius, 
mutual gravity dominates central force and vise verse. It is possible to rewrite Eqs. (|]) as 
mass-independent forms by scaling length, time, and mass as 

X = x/rh, 

t = t^o, ^^-^ 
rh = m/Msh^ = -. 

The scaled mass fh is independent of mass of particle. In our simulations, we use these scaled 
variables and solve non-dimensional equations of Eqs. @. 

A ring system is characterized by two non-dimensional parameters, optical depth r and the 
ratio rh/2rp. The dynamical optical depth is defined as 

(5) 

where x Ly and are the domain area and total number of particles in the domain. The 
optical depth can be expressed as r ~ tx/^c (e.g., Goldreich and Tremaine 1982), where t-^ and 
tc are the Keplerian time and the mean collision time. If r is large enough (typically, r ^ 1), 
clumps tend to be formed due to collision damping and self-gravity, while the clump formation 
is inhibited by differential rotation of the Keplerian motion if r < 1, as shown later. 
The other parameter rh/2rp is 

^ o.82( ^)V3(^). (6) 

2rp ^900kg/m'^ ^O^m^ ^ ' 

where Saturn's mass, Mg = 5.69 x lO^^kg, is assumed and p represents material density of a ring 
particle. Since 7r(2rp)^ and vrr^ express geometrical and characteristic gravitational cross sec- 
tions for small random velocity, this parameter regulates which mechanism dominates velocity 
changes, gravitational scatterings or direct collisions. If rh/2rp ^ 1, gravitational scatterings is 



more effective than direct collisions and vice versa (Ohtsuki 1999). For the B-ring, rh/2rp ~ 1, 
so that both gravitational scatterings and direct collisions are important. On the other hand, 
ratio rh/2rp is also related to the efficiency of the self-gravity against the tidal force of the 
central planet. Ohtsuki (1993) showed that the tidal force inhibits formation of persisting ag- 
gregates bounded by self-gravity if rh/2rp < 1. From Eq. (H), the condition rh/2rp > 1 holds 
outside A- ring. The systems with the same rh/2rp have the same dynamical condition, even 
if ao and p are different. We examine the dependence of the wake and clump formations on 
rh/2rp rather than on ao and p. 

In our simulations, we set up initial conditions as follows. The positions of particles are 
randomly distributed in the simulation region, avoiding overlapping of particles. The velocities 
except for the shear velocity of individual particles 3QoXi/2 are chosen randomly in the limited 
range from to Srh^^o, but when the initial random velocity yields Qmitiai < 2 we adopt random 
velocity large enough to be Qmitiai > 2. The density structure obtained by our simulations is 
independent of choice of random numbers used to generate these initial conditions. The initial 
velocity is immediately relaxed and the equilibrium velocity state is established within a few 
Keplerian periods. 

Equations (0) is integrated with the fourth-order Hermite scheme, which is one of predictor- 
corrector integrators and needs time derivatives of the mutual gravitational force as well as the 
mutual gravitational forces instead of past few data of position and velocity for interpolation 
(Makino and Aarseth 1992). This integrator is easily implemented to the simulations in which 
there are discontinuous phenomenon such as direct collisions and jumps due to the periodic 
boundary conditions, because the Hermite integrator does not need data of past steps. In the 
integration, most expensive part of integration is the calculation of the mutual gravitational 
forces and its time derivatives whose calculation cost is 0{N'^). To reduce computational time, 
mutual force and its time derivative are calculated by IIARP-2, which is a special purpose 
hardware for calculating gravitational force (Makino et al. 1993). "HARP" means Hermite 
AcceleratoR Pipeline and is one kind of GRAPE system ("GRAPE" means GRAvity PipE; 
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Sugimoto et al. 1990). This hardware is connected to a host computer through a communication 
interface. The host sends position, velocity, and mass of all particles to the hardware. The 
hardware calculates gravitational forces with pipelines and returns them to the host. The host 
integrates orbits utilizing the returned gravitational forces. By using this hardware, calculation 
cost is significantly reduced. 
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2.2 Calculation of mutual gravitational force and direct collision 

As stated above, the mutual gravitational force between ring particles and the direct (inelastic) 
collision play an important role in the evolution and the kinetic behavior of particle system, and 
the wake and clump formations. An equilibrium state is established by the balance between 
viscous energy gain (heating) through gravitational scatterings and direct collisions and energy 
dissipation (cooling) through inelastic collision. 

In the original Wisdom and Tremaine's method, the mutual gravity was ignored or was 
only treated as vertical force by enhancement of vertical frequency. Salo (1992b, 1995) and 
Richardson (1994) extended simulations exactly including mutual force. Richardson (1994) 
used the tree method while Salo (1995) carried out direct sum of forces exerted on a particle 
from neighbor particles within the distance -Rmax- The cut-off length -Rmax is often substantially 
less than width of simulation region. In our simulations, we adopt direct sum of forces from all 
particles with some adjustment explained below. 

Salo (1995) studied the influence of cut-off distance -Rmax and found an optimal value. But 
the larger i?max may be needed when the strong wake-like structure appears and coherent motion 
in the structure is important. We consider gravitational forces exerted on a particle from all 
other particles used in the simulation, carefully treating particles near boundary as well as in 
Salo (1995). For the effective use of HARP, we introduce the following subregion method. As 
in Fig. 1, we consider an original region (a box with thick line) and its copies. Because of the 
Keplerian shear velocity, radially inner and outer boxes have different angular velocity so that 
these boxes slide upward and downward with the shear velocity 3r2o-^a;/2. A particle at (x, y, z) 
in the central box with the velocity {vx, Vy, v^) has its images at (xinL^., y±mLy^3i}onLxt /2, z) 
with the velocity {vx, Vy ^ 3flonLx/2, Vz), where n and m are integers. When a particle in the 
central box outgoes from a certain boundary, its image at the neighbor box comes from the 
opposite boundary. We divide the simulation region into nine subregions, represented as broken 
line in Fig. 1. For each subregion, we make virtual regions with the same size of original. In 
the virtual regions, the subregion is centered. Dark and light shaded regions in Fig. 1 are 



10 



the example of a subregion and its virtual region. For particles in the subregion, we count 
gravitational forces of all particles in its virtual region. In this method, it happens that force 
from closer particle is neglected but that from more distant particle is counted. However, this 
asymmetry occurs in the in the outer part of the virtual region and does not affect the results 
because the size of the simulation is considerably larger than -Rmax Salo (1995) used. Actually 
we confirmed it performing the simulations with different size of simulation region. 

In our simulations, we adopt a hard-sphere collision model that is commonly used in previous 
simulations (Wisdom and Tremaine 1988, Salo 1991, 1992a, 1992b, 1995, Richardson 1994). A 
collision changes only impact velocity in normal direction. Normal and tangential components 
of relative velocities of colliding particles after the collision, v'^ and are described by 

V[ = Vt 

where e is the restitution coefficient of ring particles. Generally, e depends on its material 
properties as well as the impact velocity. In spite of many efforts (Bridges et al. 1984, Hatzes 
et al. 1988, Dilley 1993, Supulver et al. 1995), we have had no precise knowledge on e yet. In 
this study, we will treat e as a parameter. 

Detection of collision is one of difficulties in treating collisions. Ideally, collisions should 
be detected at the instance when colliding pairs just touch each other but it is practically 
impossible. Such problems were discussed by Richardson (1994) and Salo (1995). In this 
study, we use a simple method for detection of collision. Collision is detected as overlapping of 
colliding particles approaching each other. We search colliding pairs by using a neighbor list of 
each particle which is calculated by HARP-2 as well as mutual forces and its time derivatives. 
If a pair in the neighbor list overlaps, rij < r^^i + rpj, and is approaching, Vij ■ Vij < 0, where 
Vij and Vij are relative distance r^j = Vj — rj and velocity Vij = Vj — Vi and rp j is radius of 
particle we assume that these particles collide and velocity of colliding pair after the collision 
is calculated by Eqs. (0). If rij < r^^i + r^j but Vij ■ Vij > 0, we let the particles separate. 
However, colliding bodies can happen to overlap significantly without having been detected. In 
this case, the bodies keep sinking, which stops calculation (Richardson 1994, Salo 1995). To 
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avoid this trouble, we remove overlapping of colliding bodies by the method used by Richardson 
(1994). In his method, overlapping bodies are moved outward along the line connecting each 
center to the location where these particles are just touching. Prom conservation of center of 
mass of colliding particles, adjusted positions are given as 

x'j = Xj + Ar, ^ ' 

where the prime denotes position after correction and 

Ar = riAl^^M^Iiir,.. (9) 

We always apply this adjustment when collision is detected. 

The penetration depth of colliding bodies can be controlled by changing accuracy of timestep. 
In our simulations, we adopted shared, variable timestep because of frequent direct collisions 
and mutual gravitational forces. The formula of timestep is simply given by 

5t = min(a|^) (10) 

where a is timestep coefficient and Oj and are acceleration of particle i and its time derivative. 
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3 Results of simulations 
3.1 Comparison with Salo(1995) 

In this study, we restrict our simulations to identical particle systems in order to understand the 
basic dynamical properties. Previous simulations performed by Salo (1992b) and Richardson 
(1994) included the effect of a mass (size) distribution of ring particles so that our simulation 
results can not be compared directly with their results. On the other hand, Salo (1995) per- 
formed simulations of identical particle systems with a wider range of parameters and obtained 
statistical quantities such as radial velocity dispersion. First we perform simulations with the 
same parameters as used in Salo (1995) and confirm his results. 

Salo (1995) found the spontaneous formation of the wake structure in identical particles 
if self-gravity of ring particles is effective as well as inelastic collision, while previous studies 
without the self-gravity by Wisdom and Tremaine (1988) and Salo (1991, 1992a) did not find 
such a structure. The structure is non-axisymmetric one such as collective Julian- Toomre 
wakes (Julian and Toomre 1966). The formation of the wake structure should be related with 
gravitational instability. Salo (1995) performed simulations with various ring parameters such 
as r, e, Oq, and p (last two correspond to rh/2rp) and showed that the formation of the wake 
structure is regulated by the condition of the Toomre parameter Q value for gravitational 
instability (Eq. (|T])). Small e and large r lead to gravitational instability, because small e 
results in small velocity dispersion and large r results in large surface density. For small 
rh/2rp, formation of persistent aggregates by accretion are suppressed by the tidal force. 

Formation of spatial structure influences the motion of individual particles. Salo (1995) 
showed that as the wake structure grows, the radial velocity dispersion increases with a large 
magnitude of its fluctuation and tends to keep the relation Q ~ 2. Salo (1995) also discussed 
the equilibrium radial velocity dispersion. He showed that when rh/2rp is small, i.e., the self- 
gravity is insignificant, the velocity dispersion is dominated by direct collisions rather than the 
self-gravity to have a value of an order ~ r^Vt, but when rh/2rp is large but the surface density 
is small enough to avoid gravitational instability, gravitational encounters dominate and the 
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velocity dispersion attains the surface escape velocity v^s = y'^Gm/rp, where m is mass of a 
particle. Thus, equilibrium velocity may be altered by the parameter which characterizes a ring 
system as well as inelasticity of a ring particle. 

Our simulations show quite similar results to Salo (1995). Figure 2 is the spatial distribution 
of particles viewed face on and edge on. In Fig. 2a and 2b, the optical depth r is changed with 
fixed, velocity-independent coefficient restitution e = 0.5 and e is changed with fixed r = 0.4, 
respectively. The other parameters are the same in all simulations; rh/2rp = 0.82 corresponding 
to, for example, the situation with internal density po — 900kg/m^ and Saturncentric distance 
Qq — 1.0 X lO^m, and — Im. The simulation region is square with the width L — 112m 
{L = 68.3rh). The optical depth is controlled by changing number of particles, for examples, 
we use N — 400 for r = 0.1 and N — 2400 for r = 0.6, respectively. Fig. 2c is the results of the 
simulations with various ratio rh/2rp for r = 0.4 and e = 0.5. In these simulations, the region 
is also square and the width of the region and number of particles are L = 59.3rh {L = 80m) 
and N = 810 for rh/2rp = 0.675, L = 68.3rh {L = 112m) and N = 1600 for rh/2rp = 0.82, and 
L — 70.0rh {L — 140m) and N — 2496 for rh/2rp = 1.00, respectively. These figures correspond 
to Fig. 10b, 11c, and 9c (and Fig. 14) in Salo (1995), respectively. Figure 2 shows that the 
increase of r and/or ratio rh/2rp, and the decrease of e lead to the formation of the strong 
wake structure, as pointed out by Salo (1995). This wake-hke structure is formed after a few 
Keplerian periods from the beginning of the simulation. The wake structure changes with time. 
It is formed and dissolved successively on time scale of an order of a Keplerian period, so that 
the evolution of spatial structure affects the equilibrium velocity dispersion. 

Figure 3 shows the equilibrium velocity of (f^)^^^ and {vlY^^ as functions of r and e for 
the simulations of Fig. 2a and 2b. In Fig. 3b, we also show the results with r = 0.1. 
Each equilibrium value and error bar are the time average and the fluctuation of each velocity 
component after initial relaxation (see Fig. 4). In simulations of Fig. 2 and Fig. 3, the initial 
relaxation time is a few Keplerian periods. As shown in Fig. 4, the dispersion {vlY^^ in the 
structured case (r = 1.2) oscillates regularly with time. In this case, the bars represent the 
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amplitude of the oscillation. The dotted and dashed lines in Fig. 3a are the critical velocity 
calculated from Toomre's Q being Q — 1 and Q — 2, respectively. Table I and II list the 
equilibrium values of {v'^^^'^ scaled by r-^O-o obtained by Salo's and our simulations with various 
r and e. In Fig. 3, filled circles and open triangles are our results and Salo's, respectively. They 
are in good agreement with each other. 

Prom the comparison of Fig. 3 with Fig. 2, it is evident that the formation of the wake 
structure leads to the large increase in the equilibrium radial velocity dispersion with the large 
amplitude of its fluctuation, as pointed out by Salo (1995). On the other hand, any drastic 
change is not observed in vertical velocity dispersion {vlY^'^ even if a strong wake structure 
is formed. Figure 3a shows that (f^)^^^ increases with r, following the line of Q = 2. This 
tendency of (t'^)^^^ to follow Q ~ 2 is commonly observed in simulations with other parameters 
which allow the formation of the wake-hke structure. Figure 3b shows that in the region e ^ 0.6 
for r = 0.4 and e ^ 0.4 for r = 0.1, {vlY^"^ is reduced as decreasing e because the coUisional 
dissipation becomes strong. Thickness of ring is also reduced as clearly seen in x-z plot of 
particle distribution for e = 0.6 and 0.7 in Fig. 2b. On the contrary, in the region e ^ 0.6 for 
r = 0.4, {vl)^/"^ increases with decrease in e. In the case of r = 0.1, r is too small to cause 
gravitational instability, so that the increase of {vlY^"^ is not found. In this case, gravitational 
scatterings dominate and the equilibrium velocity is given by the surface escape velocity, as 
pointed out by Salo (1995). This tendency is also observed in Fig. 3a at small r. Note that in 
the case where rh/2rp is small, the equilibrium velocity should be dominated by direct collision 
rather than by gravitational scattering. Similar results were obtained by Ohtsuki (1999), who 
investigated the equilibrium velocity dispersion in the dilute ring system (r <^ 1) by numerical 
three-body integration of particle's orbits taking into account finite size of particles. When 
e is sufficiently large (e > 0.7 for r = 0.4, e > 0.6 for r = 0.1), coUisional damping cannot 
equilibrate with stirring due to gravitational scatterings and direct collisions, so that both 
{v^Y^"^ and {v'lY^'^ monotonically increase. These results approximately agrees with the critical 
restitution coefficient for the existence of a steady state obtained by Goldreich and Tremaine 
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(1978) and Ohtsuki (1992, 1999). 

Thus, our simulations as well as Salo (1995) show that the formation of the wake structure 
occurs in a dense ring system with self-gravitating particles and that as the spatial structure 
develops, the radial velocity dispersion increases with a large magnitude of fluctuation and takes 
Q ~ 2. Salo (1995) suggested that such increase of the velocity dispersion is due to scatterings 
by collective wakes. However, he did not give detailed argument. On the other hand, Salo 
(1995) suggested through the velocity field of particles in the wake that the velocity increase 
comes from the difference of systematic motions between adjacent wakes. Next we will analyze 
the motion of particles in detail to clarify the relation between the velocity increase and the 
wake structure. 
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3.2 Coherence of particle's motion in the wake structure 

As shown in the preceding subsection, the formation and the evolution of the spatial structure 
and the motion of particles influence each other. To clarify the relation, we analyze motion of 
particles in the wake-like structure in detail. 

First we compare the evolution of the radial velocity dispersion in the cases with and without 
wake-like structure. We perform two large A'"-body simulations for r = 0.3 and r = 1.2, 
which correspond to the non-structured and the structured cases, respectively (see spatial 
distributions in Fig. 4) . The number of particles and the width of the square simulation region 
are N = 5052 and L = 230m (L = 170.5rh, L/Acr = 16.54) for r = 0.3, and L = 115m 
(L = 85.2rh, -L/Acr = 2.07) for r = 1.2. The other parameters are fixed as e = 0.4 and 
rh/2rp — 0.675. Figure 4 shows the evolution curves of the velocity dispersion {v'^y^'^ and 
the spatial distributions in two cases. For r = 0.3, the wake structure is very weak because 
of low surface number density and (f^)^^^ attains a steady state with a small magnitude of 
fluctuation after a few Keplerian periods. On the other hand, for r — 1.2, the velocity dispersion 
increases with oscillation with a large amplitude after initial relaxation time (~ Itx)- In the 
beginning of the simulation, {v'^y^'^ is rapidly reduced by the dissipation of random motion 
energy through inelastic coUisions, and it becomes small ~ 2.5rhilo- But, after t — Itx, {viy^"^ 
begins to increase, accompanied by oscillation. This velocity increase is closely associated with 
the formation of the strong wake structure. 

Similar results are also observed in simulations with different sets of parameters in which 
wake structure is formed. To examine the effects of the boundary conditions, we perform 
simulations with various sizes of simulation domain in the structured case, fixing r = 1.2, 
e = 0.4, and rh/2rp = 0.675 (see table III). Figure 5 shows the equilibrium velocities as a 
function of number of particles N (equivalently, the size of the simulation domain). Basic 
features are the same in all simulations, although it seems that {v'^y^'^ increases shghtly with 
the increase of (simulation domain). The time evolution curves of {v'^y^'^ also show oscillation 
with a similar amplitude and a period. 
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Salo (1995) showed similar oscillation of (f^)^^^ but the oscillation is not well resolved and 
recognized as "fluctuation". The reason is that the interval of his sampling time was large 
(~ O.S^K in Fig. 2 in Salo 1995, in our simulations the interval is O-Oldtx)- We measure a period 
of the oscillation by Fourier transform methods. Figure 6 shows the results of this spectrum 
analysis in the non-structured and the structured cases. We find a strong peak at ~ O.G^k for 
T = 1.2 but not for r = 0.3. Figure 7a is the results with various computational region for 
r = 1.2 and show similar peak period. Figure 7b, 7c, and 7d are the results of simulations for 
various r, e, and rh/2rp, with the other parameters fixed. The dependence of the period on r 
and e is small. But the change of ratio rh/2rp affects the period of oscillation. Figure 7d shows 
that the width of peak and the period become larger with the increase in rh/2rp. Simulations 
in the structured cases always exhibit such an oscillatory feature and it would be related with 
the property of motion of particles. 

We analyze the concurrent evolution of the spatial and the velocity distributions in detail. 
Figure 8a and 8b show the time sequence of the spatial distribution of particles from t — 8.91iK 
to 10.19iK and the corresponding velocity distribution in the case of r = 1.2. 

As seen in Fig. 8a, temporal clumps are formed due to gravitational instability and dissolved 
due to shearing motion. The change of spatial structure looks like turbulent fluid, where vortices 
are formed and dissolved successively. Velocity distribution as well as spatial distribution is 
not steady. Figure 8b shows that the velocity distribution revolves clockwisely, holding bar-like 
conflguration. The conflguration is similar at t — 8.91, 9.55, and 10.19iK, so that the period of 
the revolution is about IMk- The oscillation in the evolution curve of {viy^^ in the structured 
case observed in Fig. 4 would come from revolution of ^^-distribution. The periods obtained 
by the Fourier analysis of the evolution curve of {v^y^'^ (Fig. 7) and in the velocity distribution 
are different by just a factor 2. This is why {v'^y^'^ is mean square of v^, i.e., the velocity 
dispersion takes maximum value twice during a period of the revolution of the configuration 
in velocity space. On the other hand, such change in the velocity distribution is not observed 
in the non-structured case. Figure 9 shows the velocity distribution for r — 0.3 at t — IAMk- 
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The distribution is independent of time. 

In the theoretical studies (e.g. Goldreich and Tremaine 1978, Araki and Tremaine 1986), 
distribution function was assumed to be Gaussian distribution. In the optically thin case 
(r <C 1), A^-body simulation showed that velocity distribution is Gaussian (Ida and Makino 
1992). The cumulative distributions of at different three time oi t — 9.55, 9.71, and 9.87iK 
in the structured case with r = 1.2 are shown in Fig. 8c and that sit t — 11.1, 12.7, and 
14.3tK in the non-structured case with r = 0.3 are shown in Fig. 9. The velocity distribution 
is always given by Gaussian distribution in the latter case while it deviates from Gaussian 
distribution in the former case. With the Kolmogonov-Smirnov test (see, e.g.. Press et al. 
1986), we calculate the probability which represents the similarity between a distribution and 
Gaussian distribution. The probability is P = 0.64 at t — II.I^k, P — 0.38 at t — 12.7iK, and 
P — 0.41 at t — 14.3iK, respectively, in the latter case. If P = 1, two distributions are the 
same but if P ^ 1 they are different. In the structured case, P = 8.18 x 10^^'^ at t = 9.55tK, 
P = 1.52 X 10-20 at t = 9.71iK, and P = 5.52 x 10"^° at t = 9.87iK, respectively. Thus 
our results show that in the non-structured case the velocity distribution of Vx is expressed by 
Gaussian distribution with high accuracy but in the structured case the distribution is far from 
Gaussian. 

How do individual particles behave in the wake-like structure? The x-component of position 

and the velocity of a Keplerian particle i relative to the local shear velocity are written as 
(Nakazawa and Ida 1988) 

Xi = aobi - aoCi cos(Qo^ - Si), 

Vx,i = aof^oCi sin(Qo^ - (^i), (H) 
'^y,r,i = -aofioCi cos(Qo^ - (^i), 

where Cj and 5j are the eccentricity and the longitude of perihelion, and ao(l + hi) is the 
semimajor axis of particle i. Vy^r,i denotes ^/-component of velocity subtracted the local shear 
velocity, Vy^r,i = Vy^i + 3QoXi/2. If motion of particles is true Keplerian, orbital elements such 
as Cj, 5j and 6j are constants and the evolution curve of x {v^) of each particle represents a 
regular oscillation in x (vx) with a period of I^k and with an amplitude of eoo (eJloOo)- Mutual 
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interactions such as direct collisions and self-gravity should cause motion of particle to deviate 
from true Keplerian orbit. Figure 8d and 8e show the time evolutions of x-components of 
position and velocity of arbitrary chosen three particles from t = IGtx to t = SOtx in the 
structured case used in Fig. 8a. On timescale of a few ^k, regular oscillation which comes from 
Keplerian epicyclic motion is seen in both x and Vx but the period of this oscillation is slightly 
larger than Itx- On longer timescale, the motion deviates from purely Keplerian orbit. We call 
such motion as "pseudo" -Keplerian motion. 

The revolution of the velocity distribution with a period of an order of a Keplerian pe- 
riod observed in Fig. 8b would reflect "pseudo"- Keplerian motion of particles in the wake. 
From Eqs. (|11|), if perturbation by other particles is negligible, trajectory of a certain parti- 
cle in velocity space is a ellipse with 2:1 axis ratio elongated to the fa;-direction expressed as 
{vx,i/cio^oGi)'^ + {vy^r,i/{^0'0^o^i)Y = 1 ^^d the particlc revolves in the clockwise direction. If 
all particles have similar phases 6.i or 6i + TT, phase space motion would look like that in Fig. 8b. 
This means that particles move coherently in the strong wake structure, holding the property 
of "pseudo" -Keplerian motion. Self-gravity of particles may force their coherent motion with 
the mechanism proposed to maintain the elliptic ring around Uranus proposed by Goldreich 
and Tremaine (1979). If the motion is true Keplerian, the oscillation period should be equal 
to li^K- But in the above structured case, the period observed in our simulations is about 
l.Stx; apparently longer than Itx- The increase in the period may be due to friction caused by 
inelastic collisions. 

Next we analyze spatial distribution of particles, comparing it with the velocity distribution. 
In Fig. 10, we mark particles in 3 local spatial regions by filled circles, triangles, and squares 
and also mark the corresponding particles in the velocity space, both in non-structured and 
structured cases. Spatial distributions used here are the results in Fig. 4. In the non-structured 
case with r = 0.3, the marked particles are scattered randomly in the velocity space, which 
means that there is no correlation of the velocity among neighbor particles, i.e., particles move 
almost randomly. But for r = 1.2, as seen in Fig. 10b, the neighbor particles in the real space 
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are also localized in the the velocity space. This indicates that neighbor particles have a similar 
velocity. Thus, in the wake, particles no longer move randomly but coherently (i.e., neighbors 
move together in the same direction at a similar velocity). 

The coherent motion of particles is gradually destroyed as the neighbor particles are sepa- 
rated from each other by shearing motion. We found that the marked particles diffuse in both 
space and velocity space on the time scale of an order of a Keplerian period. But a new coherent 
group is created from diffused distribution by the gravitational instability. Continuous creation 
and disruption of grouping occur in the wake. 

The coherent motion affects the mean velocity of particles. The bulk velocity due to the 
coherent motion increases the radial velocity dispersion. Consider the population of particles 
marked by triangles in Fig. 10. In the non-structured case, triangles are distributed randomly 
in the velocity space so that the coherent motion is weak and the width of the distribution 
represents the randomness of motion. On the other hand, in the structured case triangles are 
distributed around a non-zero velocity. This means that triangles coherently move and have 
large bulk velocity. The spread of this region represents a magnitude of the local random 
velocity. We attempt to quantitatively separate the motion of a particle into the bulk and the 
local random motions following Salo (1995). Salo (1995) mainly considered the local random 
motion but did not study the bulk motion. The bulk motion is also important for understanding 
the wake structure, because the motion of particles is largely governed by the bulk motion. We 
consider both velocities in detail. We define the bulk velocity and the local random velocity 
of particle i as the mean velocity averaged from the 10 nearest particles around particle i and 
the difference from this mean velocity (Use of more particles may inchidc particles that are not 
coherently moving while less particles may cause large statistical fluctuation). These velocities 
are written as ^ 

^b,« - tt; J2 -"j^ 

1" . 12 

lOnearest 

The (total) velocity dispersion is expressed by using the bulk and the local (random) velocity 
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dispersions as 

(^a)^ + (13) 

where a denotes components of a;, y, 2;-direction. If there is no correlation of the velocity between 
neighbor particles, i.e., particles move randomly, v^,? — so that (v^)^/^ ~ (t^f)^^^, while if there 
is strong correlation, |vb,i| ~ I^mI so that (v^)-*-/^ ~ ('^b)^''^- 

Figure 11 shows the time evolution of x-component of the bulk and the local (random) 
velocity dispersions in the non-structured and the structured cases in Fig. 4. The local velocity 
dispersion is similar in both cases, initially the velocity dispersion decreases due to the collision 
damping to attain an equilibrium value with a small magnitude of fluctuation. On the other 
hand, the bulk velocity dispersion shows a quite different feature. In the non-structured case, 
bulk velocity is smaller than local velocity. This result reconfirms that in the non-structured 
case, particles move randomly but not systematically. But in the structured case, the bulk 
velocity grows and a large amplitude oscillation starts. The amplitude is far greater than the 
local velocity dispersion. In the structured case, the bulk velocity dominates the radial velocity 
dispersion. It is clear that the increase and the oscillation in the radial velocity dispersion are 
caused by bulk motion. 

The distributions in the local and the bulk velocity spaces are also studied by the same 
method as in Fig. 10, in which neighbor particles are marked by the same symbol. Figure 
12 separates the results in Fig. 10b into the bulk and the local random velocity components. 
The distribution of the local random components shows randomness of motion as in the non- 
structured case in Fig. 10a and there is no correlation between neighbor particles in the local 
random velocity. The shape of this distribution does not change with time. The local random 
velocity dispersion represents randomness of motion. The distribution of bulk component in 
Fig. 12b shows the clockwise revolution more clearly than in Fig. 8b. 
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3.3 Characteristic scale of the wake structure in self-gravitating par- 
ticles 

We also studied the characteristic scale of the wake structure. Number of particles used in Salo 
(1995) may not be enough to obtain the scale accurately. Although Griv (1998) performed local 
A'^body simulations using = 8000 self-gravitating particles but without inelastic collision, he 
did not measure the scale. We analyze the characteristic scale of the wake through simulations 
with > 5000 up to 40000 particles. 

Simulations were done for three different optical depth r = 0.3, 0.6, and 1.2 with fixed 
e = 0.4 and rh/2rp = 0.675. For the first two simulations, the simulation region is a square 
and the width and number of particles are L = 230m {L = 170. 5rh) and N = 5052 for 
T = 0.3, L = 150m (L = 118.6rh) and = 4890 for r = 0.6. For r = 1.2, we set the size 
of the simulation region and number of particles as = 320m (L^. = 240rh), Ly = 350m 
{Ly = 260rh), and A^ = 42780. We adopted the autocorrelation analysis to determine the 
scale. The autocorrelation of the distribution function of particle position on x-y plane n{x, y) 
is defined by 

/oo /"OO 
dx I dy n{x + r,y + s)n{x,y). (14) 
-oo J ~oo 

The distribution function n{x, y) is given by Dirac's (5-function 5{x) as 



n[x. 



y) = Aj25{x-x^)6iy-yi), (15) 



where A is normalization factor and is chosen as A = ^L^Ly/N for later convenience. Substi- 
tuting Eq. ( p!5D into Eq. (0), we obtain 



Corr(r, s) = — ^^X^X^y dx 6{x + r — Xj)6{x — 



Xi 



X / dy 6{y + s - yj)6{y -yi). (16) 

The integral in terms of x is calculated as 

dx S(x + r — x,j)S(x — Xi) = \ ^ ^ 
3 10 otherwise 
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and the same calculation is done for the integral of y. This means that the right hand of 
Eq. (^) except for the coefficient represents the number of pairs of particles which satisfy the 
conditions r = Xj — Xi and s = yj — i/i, simultaneously. We take an average of Eq. (|16D in terms 
of r and s in the range [r, r + Ar] and [s, s + As] as 

]^ rr+Ar ps+As 

Corr(r, s) = ^ ^ j drj c/s' Corr(r', s'), 
= "'('■• 

where np{r, s) is number of pairs of particles which satisfy r < xj — Xi < r + Ar and s < 
Uj ~ Hi < s + As. Ar and As are given by using numbers of division of simulation region 
and A^j; as Ar = L^/N^ and As = Ly/N^ and Eq. ([I8| ) is expressed as 

np(r,s) 



This formula is the same as superposition of particle's distribution viewed from each particle. 
In the case where particles distribute homogeneously in x-y plane, i.e., in the non-structured 



case. Up ~ N"^ /NxNy so that Corr(r, s) ~ 1. 

Salo (1995) also adopted the correlation function of particle position by superposition of 
particle's distribution instead of Fourier transform. Using this analysis, Salo (1995) displayed 
the wake-like structure similar to gravitational wake in the non-colliding particle's system that 
was studied by Julian and Toomre (1966). He mainly studied the influence of the boundary 
conditions on the configuration of the wake-like structure. In the analysis of the autocorrelation 
function, we adopt the same method used by Salo (1995). Figure 13 shows the autocorrelation 
functions of particle position and the particle distribution used in this analysis. For r = 0.3, 
inhomogeneity cannot be seen in particle distribution. The autocorrelation function does not 
have any clear structure without a peak near origin which is caused by temporary gravitational 
binding pointed out by Salo (1995). On the other hand, as r increases, a wave-like structure 
gradually appears. For r = 0.6 and 1.2, several clear waves form in the autocorrelation func- 
tions. The shape of the waves slightly changes with time but the characteristic scale does not 
change. Figure 14 shows the cross sections aX y = constant for r = 0.6 and 1.2. The cross 
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sections at different y overlap each other with appropriate shifts. From these figures, we found 
that the characteristic wavelengths are \Qh/r\ ~ 25 for r = 0.6 and Aob/'^'h ~ 50 for r = 1.2, 
respectively. 

Although these waves are the results from fully developed gravitational perturbation, we 
compare our results with the linear perturbation theory. We confirm that the scale of the wake 
structure is characterized by the longest wavelength Acr that causes axisymmetric gravitational 
instability in a thin disk, as pointed out by Salo (1995). In the linear theory, the most unstable 
wavelength for a rotating and non-pressure sheet is (see, e.g., Binney and Tremaine 1987) 

Ann = (20) 

where Vg is the random velocity. When the wake structure is formed, (f^)^''^ no longer represents 
random velocity because of the development of coherent motion. We use the local velocity 
dispersion {vf^y^'^ as well as (f^)^/^ for the evaluation of Eq. (pO]). On the other hand, the 
wavelength A^ is given as (Toomre 1964) 

A„ ^ (21) 

where epicyclic frequency k, is equal to the angular velocity Qo in Saturn's ring. By using 
parameters and values of the velocity dispersions obtained by our simulations, we evaluate the 
values of these wavelengths in both cases of r = 0.6 and 1.2 and the results are listed in table 
IV. The wavelengths Aun calculated from the local velocity dispersion are quite different from 
the results of simulations in all cases. For r = 1.2, Aun calculated from (f^)^^^ is two times 
lager than observed value and the dependence of r seems to be different. On the other hand, 
Acr agrees with the observed wavelength at all r. Thus, the scale of the wake structure is 
approximately determined by A^ rather than Aun- The result of this analysis reinforces Salo 
(1992b, 1995). For the B-ring, the wavelength calculated from Eq. ( pi]) with B-ring parameter 
S = 120g/cm^ (e.g., Esposito 1993) is about 83m for r = 1.0 (we assumed Vp = Im and 
p = 900kg/m^, respectively). This value is consistent with the observational scale ~ 100m 
obtained by PPS observation (Esposito et al. 1983a, Esposito 1993). 
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3.4 Velocity change along the line Q ~ 2 

Our simulations as well as those in Salo (1995) showed Q ~ 2 when the wake structure formed. 
The Q- value is calculated by Eq. (|lD, using the equilibrium value of the radial velocity dispersion 
(v^)^/^. But we showed that the motion of particles includes the large systematic motion and 
(f^)^/^ no longer represents a magnitude of the random motion in the structured case. By using 
the local velocity dispersion, we obtain — 0.38 < 2 so that gravitational instability 

is very important (although ring particles do not form persistent clumps due to small rh/2rp 
value) . 

In Fig. lib, the simulation starts with large random motion (large local velocity disper- 
sion) to avoid gravitational instability. When the local velocity dispersion is reduced by the 
coUisional dumping to the value for gravitational instability, Q{{vf^)^^'^) ^ 2, the bulk velocity 
dispersion starts to increase. The gravitational instability makes particles move coherently and 
the coherent motion dominates the radial velocity dispersion. The oscillation pattern comes 
from the coherent "pseudo"-Keplerian bulk motion. The mean value of the oscillating velocity 
becomes constant independent of time after Six- Using the mean value of the velocity disper- 
sion, we obtain Q ~ 2. For simulations with different S we also obtain Q ^ 2. That is, the 
mean velocity is larger for larger r. 

However, we have not yet understood why Q ~ 2 is attained by the systematic motion. One 
possibility is gravitational scattering between the transient wakes, as Salo (1995) suggested. 
For the simplicity, we approximate the scattering efficiency by that between spherical clumps 
with mass A^j.S, which is typical mass of a clump formed by the self-gravitational instability. 
Gravitational perturbations pump up velocity dispersion to ~ t/j^o at one encounter in the case 
of nearly circular non-inclined orbit (e.g., Ida 1990). The Q- value in this case is (see Eq. (|l|)) 

3.36GS 3.36Gs4m, ' (22) 
~ 2.4. 

Note that the dependences of Q on particle mass, S, VLq, and Mg are canceled by Acr (see 
Eq. ([21|)) and Q takes a constant value ~ 2. When E is larger, more massive clumps forms 
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and velocity dispersion is increased by the stronger gravity. This result is consistent with the 
simulation results, however, more detail and more correct treatment is needed for the scattering 
between the wakes. 
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4 Conclusion and Discussion 



In this study, we performed local iV-body simulations including the mutual gravitational force 
between ring particles as well as the direct (inelastic) coUision with identical (up to A/" ~ 40000) 
particles. Salo (1995) showed that the spatial structures (wake and clumps) arise spontaneously 
in a dense, self-gravitating, and coUisionally dumping particle system due to gravitational 
instability. He also showed that as the wake forms, the radial velocity dispersion increases 
with the relation Q ~ 2, where Q is Toomre's non-dimensional parameter. We reconfirmed his 
results (section 3.1). Furthermore, our simulations showed a regular oscillatory pattern in the 
velocity dispersion which was not found in Salo (1995). We analyzed the motion of particles in 
the wake structure in detail and clarified that the intrinsic physics of the regular oscillation is 
coherent "pseudo"-Keplerian motion induced by self-gravitational instability (section 3.2). 

The wake-like structure is not steady but changes with time, continuously forming clumps 
due to gravitational instability and dissolving due to shear motion, which looks like the turbu- 
lent fluid. We clearly showed that in such circumstances, particles no longer move randomly 
but coherently. The coherent motion results in transient Keplerian motion with similar orbital 
elements of the particles in a clump, presumably because of mutual forcing by self-gravity. 
Such a forcing mechanism is proposed to explain the elliptic ring around Uranus (Goldreich 
and Tremaine 1979). If we separate the particle motion into the bulk and the local random 
motions following Salo (1995), in the structured system, the motion is dominated by the bulk 
motion and the local random motion is small with almost steady distribution. Owing to the 
coherent motion, the radial velocity dispersion oscillates in a period of the order of the Keple- 
rian period, but the period observed in our simulations is somewhat larger than one Keplerian 
period, which may be because of friction through inelastic collisions. 

We obtained the wavelength of the wake through autocorrelation analysis in section 3.3. The 
wavelength observed in our simulation is approximately given by the longest wavelength Apr in 
the linear theory of axisymmetric gravitational instability in a thin disk. The wavelength Acr 
is given by Acr = 47r^GE/K^ (Toomre 1966). The gravitational scattering between the clumps 

28 



with typical mass ~ EA^j. may be responsible for Q ~ 2 as shown in section 3.4, however, we 
need more detailed analysis on non-linear effects of the strong wake structure. 

The structure observed in our and Salo's simulations is non-axisymmetric one. The obser- 
vation in Saturn's A-ring produced the evidence of the asymmetric structure which may be 
caused by the wake (Dones et al. 1993). The characteristic scale of the structure observed in 
our simulations corresponds to the size of an order of ~ 100m for B-ring's parameters, which is 
consistent with the scale of the microstructure observed by PPS observation. The microstruc- 
ture may correspond to the wake structure we found. More information about structure of the 
ring will be obtained by Cassini mission. 

Our simulation results could be applicable to other optically thick particulate disks in which 
both the mutual gravitational force between particles and the coUisional dumping are important. 
One of such systems is the Uranian elliptic ring system. Uranian rings are also optically thick 
system (e.g., Esposito 1993) and a wake-like structure may arise due to gravitational instability. 
In the framework of a giant impact model of the origin of moon, it is considered that an optically 
thick debris disk formed around proto-earth after the impact and that the moon would have 
accreted in the disk (Hartmann and Davis 1975, Cameron and Ward 1976). According to 
A^-body simulations performed by Ida et al. (1997) and Kokubo et al. (1999), spiral patterns 
develop in the early stage of the simulation and it results in rapid transfer of angular momentum 
and mass. 

In the next paper, we will address angular momentum transfer enhanced by gravitational 
torque of the wake structure in detail. 
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TABLE 

Table I. Equilibrium radial velocity dispersions in simulations with parameters used by 
Salo(1995); simulations with various optical depth r for e = 0.5, rh/2rp = 0.82. 



T 


N 


L/Xcr 


Our results 


Salo(1995) 


0.1 


400 


13.5 


2.094 ± 0.084 


2.05 ±0.08 


0.2 


800 


6.7 


2.157 ±0.084 


2.16 ±0.08 


0.3 


1200 


4.5 


2.443 ±0.149 


2.46 ±0.13 


0.4 


1600 


3.4 


2.906 ±0.261 


3.01 ±0.27 


0.5 


2000 


2.7 


3.987 ±0.663 


4.05 ±0.53 


0.6 


2400 


2.2 


5.094 ± 1.137 


5.84 ± 1.42 



Table II. Simulations with various velocity-independent restitution coefficient e for r = 0.4 

and rh/2rp = 0.82. 

e Our results Salo(1995) 

0.1 4.500 ±0.589 4.39 ± 0.671 

0.2 4.578 ± 0.762 4.27 ± 0.732 

0.3 4.221 ± 0.741 3.84 ± 0.671 

0.4 4.059 ±0.589 3.48 ± 0.427 

0.5 2.906 ±0.261 3.05 ± 0.427 

0.6 2.584 ±0.168 2.61 ±0.183 

0.7 3.670 ±0.114 3.60 ±0.073 



Table III. Simulations with various number of particles for r = 1.2, e = 0.4, and 

rh/2rp = 0.675. 



N 


Lx/rh 




Lx/ Acr 


Ly/ Act 


Vx/rh^o 


err 


2460 


60 


60 


1.45 


1.45 


5.45 


1.61 


5054 


85 


85 


2.07 


2.07 


5.28 


1.49 


10696 


60 


260 


1.45 


6.27 


6.12 


1.51 


21290 


120 


260 


2.90 


6.27 


7.16 


1.38 


42780 


210 


260 


5.80 


6.27 


6.55 


0.712 
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Table IV. Characteristic wavelength for large computational simulations. 



T 




Kv/rh 










0.3 




10.31 


22.67 


15.56 


1.72 ±0.0267 


1.42 ±0.014 


0.6 


~ 25 


20.61 


19.88 


5.51 


2.27 ±0.130 


1.20 ±0.037 


1.2 


- 50 


41.22 


82.09 


3.75 


6.54 ±0.712 


1.40* 



* This value is not time- averaged. 
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FIGURE CAPTION 



FIG. 1. Schematic illustration of computational domain. Original domain is a box with 
thick line. It is surrounded by eight copied domains. The shear velocities of the outside 
and inside boxes are — 3i7o-^x/2 and +3i7o-^a;/2- Computational area is divided into nine 
subareas (broken lines) and we make a virtual area for each subregion in which the subarea 
is centered. We calculate gravitational forces in the virtual area. For example, consider 
the particle denoted by a cross in the subregion (represented by dark shaded region). This 
subregion has virtual region represented by light shaded region. 

FIG. 2a. Snapshots of the particle distribution with various optical depth (r = 0.1 to 0.6). 
Upper and lower panels are face-on and edge-on distribution. In all cases, the restitution 
coefficient e and ratio rh/2rp are e = 0.5 and rh/2rp = 0.82, respectively. The width of 
the simulation region is L = 112m {L = 68.3rh). 

FIG. 2b. The same as Fig. 2a but simulations with various restitution coefficient (e = 
0.5,0.6, and 0.7) for r = 0.4. 

FIG. 2c. The same as Fig. 2a and 2b but simulations with various rh/2rp for r = 0.4 and 
e = 0.5. The width is set as L = 59.3rh {L = 80m) for r\l2r^ = 0.675, L = 68.3rh 
(L = 112m) for rh/2rp = 0.82 , and L = 70.0rh {L = 140m) for rh/2rp = 1.00 , 
respectively. 

FIG. 3. Equilibrium values of the velocity dispersion as a function of the optical depth r 
for e = 0.5 (a) and the restitution coefficient e for r = 0.4 and r = 0.1 (b). Our results 
and Salo(1995) are represented by filled circles and open triangles, respectively. For our 
results, both (f^)^''^ and {vD^/"^ are plotted in (a) but only {vlf^ in (b). For Salo(1995), 
only {v'^y^'^ is plotted in both (a) and (b). Surface escape velocity is Ves = \J'2Gm/r^, 
where m and are mass and radius of a particle, respectively. Both dotted and broken 
lines in (a) are calculated from Toomre's Q- value (Eq. (|l])) as being Q = 1 and Q = 2. 
All velocities are scaled by t^iVLq. 

FIG. 4. Time evolution of the radial velocity dispersions for r = 0.3 (broken line) and r = 1.2 

37 



(solid line), which correspond to the non-structured and the structured cases, respectively. 
In both simulations, e and rh/2rp are e = 0.4 and rh/2rp = 0.675. The number of particles 
is N — 5052 in both cases. Typical spatial distributions of both cases are superimposed. 
For r = 0.3, edge of the distribution is trimmed away and the distribution displays the 
same area as that for r = 1.2. The width of each box is 115m (85.2rh). 

FIG. 5. Dependence of the velocity dispersion on number of particles (equivalently, area of 
the computational region). In each simulation, r = 1.2, e = 0.4, and rh/2rp = 0.675. 

FIG. 6. Spectrum analysis of a period of oscillation observed in the evolution curve of (v^)^^^ 
in the non-structured (broken line) and the structured cases (solid line). 

FIG. 7. Influence of simulation parameters on a period of oscillation observed in the evolution 
curve of (t'^)^^^. (a) Various numbers of particles but the same parameters of r = 1.2, 
e = 0.4, and rh/2rp = 0.675. (b) Various optical depth but the same parameters of 
e = 0.4, and rh/2rp = 0.675. Simulations are done with the same simulation region as 
L — 59.3rh {L — 80m) and optical depth is controlled by changing number of particles. 
Used particles in each simulation are N — 1630 for r = 0.8, N — 2038 for r = 1.0, 
and N = 2444 for r = 1.2, respectively, (c) Various restitution coefficients but the same 
parameters of r = 1.1 and rh/2rp = 0.675. Simulation region is the same as the case of 
(b), but number of particles is = 2240. (d) Various r^/2rp but the same parameters of 
r = 0.6 and e = 0.5. Width of simulation region and number of particles are L = 68.2rh 
(L = 112m) and N = 2400 for rh/2rp = 0.82 and L = 70.0rh {L = 140m) and N = 3744 
for rh/2rp = 1.00, respectively. 

FIG. 8a. Time sequence of spatial distribution of the structured case (r = 1.2) from t = 
8.91iK to t — 10.12iK- The width of region and number of particles are L — 115m 
(L = 85.2rh) and N = 5052, respectively. 

FIG. 8b. Time sequence of the corresponding velocity distribution to Fig. 8a. For each par- 
ticle, the local shear velocity {0,3floXi/2,0) is subtracted. Arrows represent the direction 
of major long axis of distribution at each time. 
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FIG. 8c. Cumulative distribution of the radial velocity Vx in Fig. 8b. Broken line with 
open circle, dashed-dotted line with open triangle, and dotted line with cross denote 
distributions oX. t — 9.55tK, t — 9.71tK, and t — 9.87tK, respectively. Sohd hne denotes 
cumulative Gaussian distribution P{s) = (1 + erf(s/\/2))/2, where erf(s) is the error 
function. The distributions are scaled by the velocity dispersions at each time. The 
instantaneous velocity dispersion is (f^ ) ins /rhf^o = 6.100 t — 9.55tK, ('f^x)ins/^h^^o = 
7.115 at t = 9.71%, and {vDH^^ / tyiVLq = 4.033 at t = 9.87tK, respectively. The time 
averaged velocity dispersion is {vDHl/rh^o — 5.862. The cumulative number is scaled by 
the total number of particles. 

Fig. 8d. Time evolutions of x-component of position (left) and velocity (right) for certain 
three particles in the structured case (r = 1.2). Solid, broken, and line-dotted lines 
illustrate the evolution curves of particles. Dotted lines denote positions of boundary of 
x-direction, the range of simulation domain is [— 42.62rh, 42.62rh]. In this figure, curves are 
connected continuously by taking into account boundary conditions. But in fact, evolution 
curve for a particle over a boundary discontinues there and the curve corresponding to 
incoming particle starts from the opposite side. 

Fig. 8e. Time evolutions of for particles used in Fig. 8d. Each origin moves with 
appropriate shift and each dotted line denotes the line oi — 0. 

FIG. 9. Velocity distributions of the non-structured case (r = 0.3). (a) Particle distribution 
in velocity space at t — 14.3^k- (b) Cumulative distribution of radial velocity scaled by 
{viy^^. The velocity dispersion is {viy^^/rh^o — 1.720. Broken hne with open circle, 
dashed-dotted line with open triangle, and dotted line with cross denote distributions 
at t — ll.ltK, t — 12.7iK, and t — H.SIk, respectively. Solid line denotes cumulative 
Gaussian distribution. 

FIG. 10. The relation of marked particles in the real and the velocity space in the non- 
structured case (r — 0.3) (a) and structured case (r — 1.2) (b). Particles within certain 
three regions with the radius 5rh which are located at (— 34rh,llrh), (10rh,30rh), and 
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(23rh,3rh) in the real space are marked by symbols as filled triangle, circle, and square, 

respectively. The other particles are represented by small dots for convenience. In the 

non-structured case, edge of the spatial distribution is also trimmed. 
FIG. 11. Evolution curves of the bulk and local velocity dispersion in the non-structured 

case (a) and the structured case (b). Dotted line , dashed hue with filled triangle, and 

solid hue denote local, bulk velocity dispersion, and (t"^)^^^, respectively. 
FIG. 12. The Same as Fig. 10b but for particle distributions in the local and bulk velocity 

space, (a) Distribution for local velocity, (b) Distribution for bulk velocity. The local 

shear velocity of each particle is also eliminated. 
FIG. 13. Autocorrelation functions of particle positions in the case of r = 0.3 (a), 0.6 (b), 

and 1.2 (c). In all simulations, e = 0.4 and rh/2rp = 0.675 are set up. We also showed 

the original snapshots used in this analysis. 
FIG. 14. Cross sections of autocorrelation functions for r = 0.6 (a) and r = 1.2 (b) in 

Fig. 13. Each cross section is shifted in the x-direction for peak positions to coincide with 

each other. 
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